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ABSTRACT 

During  the  annual  Gray  Whale  migration  from  the 
Aleutians  to  Baja  California,  the  mammals  travel  in  coastal 
waters,  thereby  presenting  an  opportunity  for  the  study  of 
their  sound  spectral  and  source  levels  and  vocabulary.   How- 
ever, such  measurements  are  distorted  by  surface  and  bottom 
reverberation,   Using  the  theory  of  rough  surface  scattering, 
knowledge  of  the  bottom  impedance,  and  correlation  techniques, 
it  is  possible  to  decompose  the  shallow  water  reverberation 
into  the  contributions  from  different  paths .   From  this ,  the 
range,  depth  and  the  deverberated  spectral  source  levels  of 
the  sounds  of  the  mammal  can  be  determined  by  use  of  only  one 
hydrophone  rather  than  the  conventional  three  or  four.   The 
theory,  deverberation  programming,  and  experimental  results 
are  presented  for  a  model  of  the  whale's  pulsed  radiation  in 
a  laboratory  model  coastal  environment. 
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I,   INTRODUCTION 

Once  yearly,  the  Gray  Whale,  Eschrichtius  glaucus , 
migrates  southward  from  the  Aleutians  and  passes  very  close 
to  the  California  coast  in  shallow  water.   During  this 
migration,  the  sounds  in  the  water  close  to  the  whales  can 
easily  be  recorded;  but,  they  may  not  be  the  true  sounds 
produced  by  the  whale.   The  whale  normally  produces  an 
intermittent,  short  duration  signal  which  in  shallow  water 
is  received  at  the  hydrophone  via  direct,  surface  reflected 
and  bottom  reflected  paths.   Since  the  path  lengths  are 
different,  the  signals  arrive  at  the  hydrophone  at  different 
times  and  they  interfere.   To  obtain  the  true  sounds  produced 
by  the  whales,  this  interference,  which  is  called  reverbera- 
tion, must  be  eliminated. 

The  purpose  of  this  thesis  is  to  model  in  the  laboratory 
the  whale  sounds  in  the  shallow  water  and  to  develop  a 
technique  to  determine  the  range,  eliminate  the  surface  and 
bottom  reverberation  and  calculate  the  spectral  source  levels 
as  a  function  of  time. 


II.   THEORY 

In  a  reverberant  environment,  the  original  signal  can 
only  be  realized  if  the  reverberation  is  removed.   The 
method  used  to  remove  the  reverberation,  which  is  here 
called  "deverberation,"  assumes  that  the  whale  is  a  point 
source,  the  geometrical  spreading  is  spherical,  the  water  is 
isovelocity  and  the  water  depth  is  constant. 

Before  the  deverberation  technique  can  be  applied,  the 
direct,  surface  reflected  and  bottom  reflected  path  dis- 
tances must  be  known.   Normally,  this  information  is  obtained 
by  knowing  the  source  position  and  calculating  the  path 
distances  from  the  known  geometry.   Determining  the  hori- 
zontal source  position  is  difficult  and  requires  at  least 
two  directional  or  three  omnidirectional  hydrophones,  accu- 
rately fixed  with  respect  to  each  other  at  all  times.   To 
determine  the  depth  requires  at  least  one  additional  hydro- 
phone at  a  different  depth.   Generally,  three  or  four 
hydrophones  are  deployed  CRefs.  1  and  21.   With  each  addi- 
tional hydrophone,  the  complexity  of  the  system  is  increased 
since  each  hydrophone  requires  its  own  amplifying  and  record- 
ing system.   In  shallow  water,  however,  taking  advantage  of 
the  surface  and  bottom  reflections,  one  can  use  a  single 
hydrophone  and  gather  all  the  information  necessary  for  the 
application  of  the  deverberation  technique. 

Consider  the  direct,  surface  scattered,  and  bottom 
scattered  sounds  received  at  only  one  hydrophone  which 
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is  deployed  in  shallow  water  as  depicted  in  figure  1.   It 
will  be  shown  that  when  the  differences  between  the  arrival 
times  for  the  three  paths  are  known  the  three  path  distances 
can  be  calculated.   Using  the  surface  specularly  scattered 
path,  R  ,  it  is  seen  that 
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Substituting  and  rearranging  terms  gives 
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Using  x  ,  the  time  difference  between  the  direct  path 
arrival  and  the  surface  reflected  path  arrival  and  C,  the 
mean  speed  of  sound,  gives 

C  Ts  -  R's  +  R"5  -  R 

and   therefore 
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Similarly  using  the  bottom  specularly  scattered  path 
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where 

Substituting  and  rearranging 

W  KB   (22  -H  -D)   lK  J 

Using  Tg ,  the  time  difference  between  the  direct  path 
arrival  and  the  bottom  reflected  path  arrival  gives 

and,  therefore, 

(6)  ^TB  =  \tt*+*U2+HD-*U  "^D^   -  R 

Solving  equations  (3)  and  (6)  simultaneously  for  R,  the 
range  of  the  source  from  the  receiver  and  D,  the  depth  of 
the  source,  produces 

°        '  ^[HtB+T5U-H)] 

(8)  R  "     2CTS 

The  equation  for  the  range  is  left  as  a  function  of  the  water 
depth  to  facilitate  its  being  programmed.   Now  that  the  range 
and  depth  are  known,  the  other  two  path  distances  can  be 
calculated.   From  equations  (1)  and  (2) ,  the  surface  reflected 
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path  distance  is 


X 


(9)  fts  =  (RZ  4  ^HD^)Z 

and  from  equations  (4)  and  (5) ,  the  bottom  reflected  path 
distance  is 

(10)  rb   =  \^2-*  ^?.z-v  WD  -  ZW-  £D)]2 

Therefore,  when  the  time  arrival  differences  for  the  differ- 
ent paths  are  known,  equations  (8) ,  (9) ,  and  (10)  can  be 
used  to  determine  the  first  three  path  distances.   The  paths 
for  multiple  reflections  can  be  calculated  directly  from  the 
known  geometry,  assuming  specular  scatter. 

For  a  transient  signal,  determination  of  the  differential 
arrival  times  t„  and  tr  can  be  realized  by  the  use  of  an 
autocorrelation  technique.   The  correlation  function  can  be 
defined  as  [Ref ,  3] 

(ii)  R  (/t)  --  E^[HtVu]lvUn)-u.]} 

where  v  is  the  time  average,  u  is  the  mean,  and  E  is  the 
expected  value  of  the  received  signal.   This  process  is 
programmed  using  a  digital  summation 

(i2)RCr)=  -^  ,S(  Kto-a^vurti}  -  u] 

with  n  representing  the  total  number  of  samples  in  the  record 
Performing  the  autocorrelation  on  the  reverberant  signal  at 
the  one  hydrophone  gives  peaks  at  delay  times  corresponding 
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to  zero  delay  time,  and  the  arrival  delay  times  from  the 
reflected  signals .   The  peaks  are  realized  only  when  the 
direct  signal  is  delayed  enough  to  correlate  with  the 
reflected  signals.   The  computer  program  called  AUTOPEAK 
performs  the  autocorrelation  and  then  applies  an  envelope 
detection  to  determine  the  delay  times  for  the  peaks .   These 
delay  times  are  then  used  in  equations  (7) ,  (8) ,  (9) ,  and 
(10)  to  determine  the  range,  depth,  and  reflected  path  dis- 
tances.  Thereby,  all  the  geometrical  information  necessary 
for  the  deverberation  technique  has  been  obtained. 

Computer  programs  have  been  developed  for  deverberation 
in  either  the  frequency  or  time  domain. 

The  computer  program  designed  to  eliminate  the  reverbera- 
tions in  the  frequency  domain  is  called  DEVERB .  For  the  time 
being,  assume  only  one  frequency  component,  whose  amplitude 
and  frequency  are  functions  of  time.  For  a  signal  with  time- 
varying  frequency  and  amplitude,  we  can  describe  the  pressure 
at  the  receiver,  due  only  to  the  direct  path  signal  as  [Ref . 
4] 

(13)  P0(O=    cU)  ^ 

Then,  taking  into  account  spherical  divergence,  the  spatial 
phase  shift,  and  specular  scattering  from  a  Gaussian  rough 
surface,  the  pressure  at  the  hydrophone  due  to  the  surface 
scattered  signal  can  be  written  as 

(u)  psw) = p9  (.t--o  --  %  e  2  e  U) 
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and  the  pressure  due  to  the  once  bottom  reflected  signal  is 

G^  and  e~^'   are  the  frequency-dependent  pressure  amplitude 
reflection  coefficients,  and  "S*  and  tt  are  the  phase  shifts 
due  to  the  bottom  and  surface  reflections,  respectively. 
The  surface  reflection  coefficient  depends  on  the  roughness 
of  the  surface.   The  exponent  for  specular  scattering  is 
[Ref.  5] 

(16)  s  *  =    1y£   C.D5  S5 

where  CT  is  the  R.M.S.  wave  height,  A  is  the  wavelength  of 
the  sound  signal  and  Go  is  the  angle  of  incidence.   Equations 
(14)  and  (15)  are  derived  from  previously  received  direct 
path  pressures,  corrected  for  path  differences  and  phase 
shifts  and  represent  the  pressures  due  to  the  single  reflected 
paths.   The  coherent  sum  of  equations  (13),  (14)  and  (15)  is 
the  pressure  sensed  by  the  receiver  when  these  three  compo- 
nents are  present. 

The  signal  processing  is  done  digitally.   Therefore,  the 
continuous  time  dependence  is  replaced,  whenever  it  occurs 
in  the  previous  equations,  by  digital  block  numbers,  indi- 
cated by  the  subscript  index  "k . "   Each  block  contains 
enough  data  samples  of  the  incoming  signal  to  give  the  desired 
spectral  frequency  resolution  during  a  block  duration  which 
is  small  compared  to  the  total  duration  of  the  time-varying 
signal.   The  frequency  change  of  the  signal  within  a  block 
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duration  is  assumed  to  be  much  smaller  than  the  frequency 
resolution  of  the  digital  processing.   The  block  duration  is 

T    NUMBER  OF  SAMPLES  IN  THE  BLOCK 
SAMPLING  FREQUENCY 

The  time  the  block  ends  is  related  to  the  continuous  time  t 
by 

(17)   t  =  TK;  K  =  1,  2,  3  • • • 

The  index  :'i"  is  used  for  the  spectral  frequency  compo- 
nent of  the  complex  wave  being  analyzed. 

The  equation  for  the  source  pressure  at  unit  distance 
using  the  frequency  deverberation  correction  is 


-[ 


■[t«.Ve'K*'*"**,Vo 


'SB 

-  ETC.} 

where  l(K-N) ,  l(K-M)  and  l(K-L)  are  unity  factors  with  values 
l(K-N)  =  1   K  >_  N 

=  0   otherwise 
l(K-M)  =1   K  >  M 

=  0    otherwise 
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l(K-L)  =1    K  >  L 

=  0   otherwise 

t""L| 

The  pressure  amplitude  of  the  i    frequency  component  in 
block  K  of  the  receiver  reverberant  signal  is  represented 
by  C^.  .  and  its  phase  represented  by  <\>„    .  .      The  deverberated 
pressure  amplitude  is  represented  by  D„  .  and  its  phase  by 

IN.  ,  J. 

^  r,  •  .   The  first  term  on  the  right  hand  side  of  the  equation 
represents  the  received  signal,   The  second  term  represents 
the  correction  due  to  a  single  specular  scatter  from  the 
surface;  the  third  represents  a  single  bottom  reflection 
correction,  and  the  fourth  represents  the  correction  for  a 
path  which  includes  one  surface  and  one  bottom  reflection. 
The  equation  can  be  expanded  to  include  other  multiple  reflec- 
tions . 

The  block  indices  are  determined  by 


(19)  M  =  K  -  ^A 
T 


(20)  N  =  K  -  -2l 

T 


(2i)  l  ■  k  -  ~ 


^  SB 


T 

The  output  of  the  frequency  deverberation  program  is  a 
series  of  spectra  of  the  Consecutive  blocks  and  its  Fourier 
transform  is  a  time  plot  of  the  deverberated  signal. 

The  above  procedure  takes  the  signal  from  the  time 
domain  (the  time  series  after  A/D  conversion)  by  Fourier 
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transform  to  the  frequency  domain,  where  the  known  frequency 
dependent  reflection  coefficients  are  easily  applied,  and 
then  back  to  the  time  domain  to  verify  the  effectiveness  of 
the  process. 

When  the  reflection  coefficients  can  be  assumed  to  be 
frequency-independent,  a  simple  point-by-point  deverberation 
procedure  can  be  applied  in  the  time  domain.   The  applicable 
temporal  deverberation  equation  is 

(22)T>lt.cK*<e"*>^DH-«%B^*<tT><R|w^*€TC- 

C^-  represents  the  pressure  amplitude  for  the  K   sample.   The 
other  terms  are  similar  to  those  in  equation  (18) .   For  low 
roughness  surfaces  g  <  1  and  the  use  of   e  °  over  the 

appropriate  frequencies  will  be  a  good  approximation  which 
is  essentially  independent  of  frequency.   One  advantage  of 
the  temporal  deverberation  technique  is  the  relative  freedom 
from  restrictions  of  block  size;  the  block  size  is  determined 
only  by  the  desired  frequency  resolution  and  the  rate  of 
change  of  frequency  of  the  transient  source. 
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III.  PROCEDURE 

In  order  to  model  the  Gray  Whale's  environment  in 
Monterey  Bay,  a  three  meter  cube  "anechoic"  water-filled 
tank  was  used.   An  artifical  bottom  made  of  hard  rubber  with 
an  experimentally  determined  pc  of  approximately  2.4  x  10 
mks  rayls  was  positioned  one  meter  below  the  water  surface. 
This  type  of  bottom  was  chosen  for  its  specific  acoustic 
impedance  since  the  bottom  at  the  listening  area  in  Monterey 
Bay  has  a  pc  of  approximately  3  x  10   mks  rayls.   The  depth 
of  the  bottom  and  relative  placement  of  the  source  and 
receiver  were  determined  in  order  to  obtain  realistic  modeled 
delay  times  between  the  reflected  signals  and  the  direct 
signal.   A  1.8  cm  diameter  spherical  source  was  used  because 
of  its  small  size  and  its  ability  to  transmit  a  signal  with 
minimum  distortion;  but  it  also  limited  the  minimum  frequency 
to  about  10  kHz,   An  FM  up-sweep  of  varying  widths  and  sweep 
rates  was  used  to  model  one  of  the  sounds  produced  by  the 
Gray  Whale. 

The  equipment  was  connected  as  in  figure  2  with  the  mas- 
ter unit  pulse  generator  being  used  to  synchronize  the  samp- 
ling frequency  oscillator  which  determined  the  start  and  stop 
frequencies,  sweep  rate  and  pulse  width.   The  pulse  repeti- 
tion rate  was  determined  by  the  master  unit  pulse  generator. 
The  signal  was  then  amplified  and  sent  to  the  source.   The 
reverberent  signal  was  received  by  an  LC-10  hydrophone,  and 
then  amplified  again  and  sent  to  the  analog-to-digital 
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converter.   The  A/D  converter  was  triggered  through  the 
delay  generator  via  the  unit  pulse  generator  and  the  sampling 
frequency  oscillator.   A  sampling  frequency  of  320  kHz  was 
used,   The  delay  generator  delays  the  trigger  by  the  time 
required  for  the  signal  to  be  transmitted  through  the  water 
and  then  this  delayed  pulse  triggers  the  unit  pulse  generator 
The  unit  pulse  generator  was  used  to  determine  the  total  on- 
time  of  the  sampling  frequency  oscillator.   The  oscillator 
determined  the  samples  per  second  taken  by  the  converter 
during  the  oscillator's  on-time.   This  complex  equipment  set- 
up was  designed  to  allow  the  A/D  converter  to  sample  only 
during  the  time  when  the  direct  and  reflected  signals  were 
present,  thereby  reducing  the  amount  of  computer  time  and 
storage  required  to  process  the  data.   Reflections  from 
other  surfaces  in  the  tank  were  eliminated  wherever  possible 
in  the  sampling  time  by  varying  the  pulse  repetition  rate 
of  the  master  pulse  generator  or  by  varying  the  source  and 
receiver  placement.   After  the  analog  signal  was  changed  by 
the  converter  to  digital  form,  it  was  stored  on  cassette 
memory  and  then  analyzed  using  the  programs  AUTOPEAK  and 
DEVERB . 
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IV.   DATA  PROCESSING  AND  RESULTS 

The  autocorrelation  plots  of  reverberation  for  two 
different  sweep  widths  are  seen  in  figure  3  with  the  scale 
factor  for  the  delay  time  equal  to  3.125  y  seconds.   The 
plot  of  the  90  kHz  sweep  width  shows  a  steeper  slope  of  the 
envelope,  thus  a  more  clearly  defined  peak  than  the  one  with 
only  a  10  kHz  sweep  width.   This  indicates,  as  expected,  that 
as  the  difference  between  the  upper  and  lower  frequencies 
decreases,  the  correlation  peaks  become  harder  to  determine. 
In  the  limit  of  only  one  frequency  being  present,  there  would 
be  no  peaks  in  the  envelope.   Figure  4  shows  the  range  error 
(computer  determination  compared  to  direct  measurement)  versus 
the  ratio  of  the  upper  frequency  to  the  lower  one.   The  graph 
indicates  that  for  a  ratio  above  1.2  to  1  the  frequency  sweep 
of  the  signal  is  sufficient  to  get  accurate  time  differences 
for  the  reflections  and  thereby  to  determine  the  source  range 
and  depth  from  the  autocorrelation  processes . 

The  frequency  deverberation  program  is  designed  to  give 
a  true  spectrum  of  the  signal  by  eliminating  frequency  depen- 
dent reverberations.   Since  the  signal  is  time-varying,  a 
small  block  of  time  is  desirable  to  keep  the  change  in  the 
signal  to  a  minimum  during  the  block.   The  limiting  factor 
to  the  minimum  block  size  is  the  desired  frequency  resolu- 
tion.  For  a  spectrum  to  accurately  represent  an  instant  of 
time  rather  than  being  a  spectrum  averaged  over  a  length  of 
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FIGURE    3. 
AUTOCORRELATION   OF   DIRECT   AND   SURFACE   REFLECTED   CHIRP 
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time,  the  frequency  change  during  the  block  time  should  be 
much  less  than  the  frequency  resolution.   In  addition,  where 
possible,  the  block  duration  should  be  equal  to,  or  a  sub- 
multiple  of  | T   -  T,  |  in  order  to  permit  equation  (18)  to 
be  applied, 

From  figure  4,  it  is  known  that  depth  and  range  can  be 
determined  only  when  the  frequency  ratio  of  the  chirp  is 
greater  than  1,2  to  1,   Figure  5  shows  the  reverberant  and 
deverberated  signals  for  a  chirp  from  46  kHz  to  54  kHz  with 
a  ratio  of  1,17  to  1.   The  reverberant  signal,  top  of  figure 
5,  was  divided  into  blocks,  equal  to  tb  -  ig,  and  then  trans- 
formed back  into  the  time  domain  which  is  shown  at  the  bottom 
of  the  figure,   Some  improvement  can  be  seen  but  the  expected 
slowly  increasing  frequency  at  the  approximately  constant 
amplitude  has  not  been  realized.   It  is  believed  that  the 
inadequate  deverberation  is  due  to  the  fact  that  the  frequency 
sweep  during  a  block  is  approximately  four- tenths  of  the 
frequency  resolution.   A  slower  sweep  rate  or  a  larger  block 
duration  would  have  cured  this  problem.   However,  a  slower 
sweep  rate  for  the  model  would  have  decreased  the  accuracy  of 
the  range  determined  by  the  autocorrelation;  and  a  large  block 
duration  was  precluded  by  the  geometry  which  determined 
|  x.g  -  To  |  .   The  temporal  deverberation  technique  which  is 
presented  next  did  not  suffer  from  these  limitations . 

The  result  of  applying  the  temporal  deverberation  program 
to  a  signal  with  a  sweep  width  from  5  kHz  to  95  kHz  can  be 
seen  in  figure  6  with  the  reverberant  signal  on  the  top  and 
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FIGURE  5. 
FREQUENCY  DEVERBERATION 
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TEMPORAL   DEVERBERATION 
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the  deverberant  signal  below,   Since  the  source  was  resonant 
at  64  kHz,  the  FM  sweep  grows  in  amplitude  to  64  kHz  and 
then  decreases.   The  deverberated  signal  shows  the  frequency 
and  amplitude  modulation  cleanly  until  the  end  of  the  direct 
signal,   The  reverberation  after  that  time  is  due  to  scatter- 
ing from  the  side  walls  of  the  tank.   Possible  reverberation 
due  to  the  fourth  and  later  terms  on  the  right  hand  side 
of  equation  (22)  were  excluded  by  limiting  the  duration  of 
sampling  by  the  computer. 
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V.   CONCLUSIONS 

The  work  described  in  this  thesis  demonstrates  the 
feasibility  of  obtaining  a  non-reverberant  spectrum  of  a 
transient  source  in  a  reverberant  environment.   The  technique 
includes  calculation  of  the  autocorrelation  of  the  received 
signal  to  determine  range  and  depth  of  the  source  and  com- 
puter processing  to  correct  for  the  surface  and  bottom  reflec- 
tions , 

The  autocorrelation  function  provides  an  accurate  method 
for  obtaining  the  range  and  depth  of  a  source  of  transient 
sound  in  shallow  water.   The  correlation  technique  can  be 
performed  for  a  chirp  sound  with  ratio  of  upper  to  lower 
frequency  of  greater  than  1.2  to  1.   At  least  two  reflections 
are  required  to  obtain  the  depth  and  range  of  the  source  with 
respect  to  the  receiver. 

Frequency  and  time  deverberation  programs  which  use  the 
position  data  from  the  autocorrelation  have  been  developed 
to  eliminate  the  reverberations  and,  thereby,  to  obtain 
corrected  spectra  or  corrected  time  plots .   Either  technique 
can  be  used;  however,  because  the  output  of  the  computer  is 
a  time  series,  it  is  natural  to  apply  temporal  deverberation. 
This  becomes  very  simple  if  the  surface  and  bottom  reflection 
coefficients  are  independent  of  frequency. 
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COMPUTER  PROGRAM  AUTOPEAK 

1  REM     DATE  OF  LAST  CORRECTION:   1C/31/77 

2  REN 

5  REN  **4##*j^5M«***#****#*#**4******#*#**5M«#****#* 

10  REM  *  AUTOPEAK   —   9/29/77  * 

15  REM  *  JEANIE  SAVAGE,  PROGRAMMER  * 

20  REM  *  SPECIFICATIONS  BY  RICK  eCSTIAN  * 

25  REM  *  THIS  PROGRAM  PERFORMS  AN  AUTO-  * 

30  REM  *       CORRELATION  ON  SIGNAL  DATA  AND  * 

35  REM  *       THEN  PICKS  OUT  THE  TfcG  PEAK  * 

40  REN  *       VALUES.  * 

45  REM      *#*#***W***#***#######*##******##**#**#*'isaic 

50  REM 

55  REN 
100  DIN  V$(48)  ,V1$(6) ,Z$( 1) 
105  DIN  A(IOOO) ,6(1000), Y( 1000) 

107  DIM  Zl$( 1) ,Z2$( 1) ,S(1) ,D(1), Z4$( 1) 

108  DIM  P(500,l)  ,Z3$( 1  )  ,M(50,2) ,R(1,2) 
110  REM 

112  Z4$="N" 

115  Z$="Y" 

115  REN  ***DRIVER  ROUTINE*** 

120  REN  PERFORM  I NI TI AL I ZAT I CN  PROCEDURE 

125  GOSUB  500  , 

127  IF  Z4$="Y"  'GCTO  140 

13C  REM  READ  DATA  FROM  TAPE 

135  GOSUB  1000 

140  REN  BUILD  CTHEP  ARRAY 

145  GOSUB  1300 

15C  REN  PERFORM  AUTO-CORRELAT ICN 

155  GOSUB  1500 

160  REM  DETERMINE  INTERMEDIATE  PEAK  VALUES 

165  GOSUB  2000 

17C  IF  Z3$="N"  GOTO  185 

175  REM  PRINT  INTERMEDIATE  VALUES 

180  GOSUB  2500 

185  REN  DETERMINE  TWO  PEAK  VALLES 

190  GOSUB  3000 

192  IF  Z8=l  GOTO  150 

205  REN  PRINT  TIME  DIFFERENCES 

210  GOSUB  4000 

215  REM  CALCULATE  SOURCE  DEPTH,  S-R  RANGE 

220  GOSUB  4500 

225  REN  PRINT  DEPTH  AND  RANGE 

23C  GOSUB  5000 

235  REM  CALCULATE  REFLECTION  PATH  DISTANCES 

240  GOSUB  5500 

245  REM  PRINT  DISTANCES 

250  GOSUB  6000 

255  PRINT  "ARE  YCU  FINISHED?  (Y  OR  N)--" 

26C  INPUT  Z$ 

265  IF  Z$="Y"  THEN  STOP 

267  PRINT  "SAME  DATA?  (Y  OR  N) — " 

266  INPUT  Z4$ 

270  PRINT  "SAME  PARAMETERS?  (Y  OR  N) — " 

275  INPUT  Z$ 

280  IF  Z$="N"  GOTO  115 

285  Z2$="Y" 

29C  GOTO  127 

295  REM 

500  REM  ***INITIALIZATION  PRGC ECURE*** 

503  IF  Z$="N"  GOTO  525 

505  PRINT  "AUTOPEAK" 

51C  PRINT 

512  PRINT  "SAMPLING  FREQUENCY  MUST  BE  FOUR  TIMES  THE" 

513  PRINT  "   GREATEST  FREQUENCY  OF  INTEREST." 

514  PRINT 

515  PRINT  "NUMBER  OF  PCINTS  PER  BLOCK--1' 

520  INPUT  N2  T  „„,  OATi 

525  PRINT  "MINIMUM  TIME  DIFFERENCE  BETWEEN  DIRttCT  PATH 
527  PRINT  "   AND  FIRST  REFLECTED  PATH  RECEPTION  (IN" 

526  PRINT  PRINT  «   SAMPLES) — " 
530  INPUT  D8 


29 


535  PRINT  "NUMBER  CF  PCINTS  TO  BE  USEC  FROM  EACH  PLCCK" 
540  PRINT  "    (NUMBER  CF  PCINTS  PLUS  DELAY  MUST  EE  LESS" 
545  PRINT  "    THAN  NUMBER  CF  POINTS  IN  BLOCK)  — " 
55C  INPUT  Nl 

555  PRINT  "MAXIMUM  LAG — " 
560  INPUT  LI 

562  PRINT   "PRINT  PEAK  VALUES?  (Y  OR  M ~ " 
562  INPUT  Z3$ 

565  IF  Z$="N"  THEN  PRINT  "CONTINUING  WITH  CALCULATIONS" 
5  7C  IF  Z$="N"  THEN  RETURN 

575  PRINT  "IS  THIS  THE  1ST  BLOCK  OF  A  MLLTIPLE  RUN?  (Y/N) — 
580  INPUT  Z2$ 

585  PRINT  "SAMPLING  FREQUENCY — " 
590  INPUT  SI 

595  PRINT  "SPEED  OF  SOUND  (M/SEC) — " 
60C  INPUT  C 

605  PRINT  "DEPTH  OF  WATER  (MJ  — " 
61C  INPUT  H 

615  PRINT  "DEPTH  OF  RECEIVER  (M)~" 
62C  INPUT  R 
62  5  RETURN 
630  REM 
1000  REM  ***ROUTINE  TO  READ  DATA  FROM  TAPE*** 

1004  PRINT  "SWITCH  TO  HIGH  SPEED." 

1005  IF  Z2$="Y"  THEN  INPUT  ON  (2)  V$ 
1010  K1=0 

1015  FOR  1=1  TO  (N2/8) 

1020  INPUT  ON  (2)  V$ 

1025  FOR  J=l  TO  48  STEP  6 

1030  Vl$="       " 

1032  Vl$(l,6)="»« 

1035  FOR  K=0  TO  5 

1040  IF  V$( J+K, J+K)="  "  GOTO  1050 

1045  V1$=V1$+V$( J+K, J+K) 

105C  NEXT  K 

1055  A(KI)=VAL(V1«) 

1060  K1=K1+1 

1065  NEXT  J 

107C  NEXT  I 

103C  PRINT  "CONTINUE?" 

1085  INPUT  11$ 

1090  PRINT  "CONTINUING  WITH  CALCULAT I CNS  ." 

1095  RETURN 

1100  REN 

1300  REM  ***BUILD  CTHER  ARRAY*** 

1305  FOR  1=0  TO  N 1-1 

1310  B(I)=A( I+D8) 

1315  NEXT  I 

1320  RETURN 

1325  REM 

1500  REM         ***CROSS-CGRRELATION  PCUTINE*** 

1505  A8=0 

151C  B8=C 

1515  FOR  1=0  TO  Nl-1 

1520  A8  =  A8+A(  I) 

1525  B8=B8+B(  I) 

153C  NEXT  I 

1535  AS=A3/N1 

154C  B8=B8/N1 

1545  FOR  1=0  TO  Ll 

1550  N9=N1-I 

1 5  5  c  S  8=  C 

1560  FOR  J=0  TO  N9-I 

1565  I9=J+I 

1570  S8=S8+( A(J)-A8)*(B(I9)-B8) 

1575  NEXT  J 

1580  Y(I)=S8/N9 

1585  NEXT  I 

159C  RETURN 

1 5  9  c  REM 

2000  REM  ***DETERMINE  INTERMEDIATE  PEAK  VALUES*** 

2005  IF  Y(1)>Y(0)  THEN  E=l 
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2010 
2015 
202C 
2025 
2030 
2035 
204C 
2045 
205C 
2055 
206C 
2065 
2070 
2075 
208C 
2085 
209C 
2095 
2500 
2505 
251C 
2515 
2520 
2525 
253C 
2535 
2540 
2545 
2550 
2555 
2560 
2565 
30CC 
3003 
3005 
301C 
3015 
302C 
3025 

3030 
3035 
3040 

3045 
305C 
3C55 
3060 
3065 
3070 
3075 
3076 
3077 
3078 
3079 
308C 
3081 
3062 
3085 
3090 
3095 
3100 
3105 
3110 
3115 
3120 
3125 
3130 

3137 
314C 
3145 


IF    Y(l)< 

P9=-l 
FOR    1=2 
IF    E=(-l 
IF    Y(  n> 
P9=P9+1 
P (P9,0)= 
P(P9,1)= 
E  =  -£ 

GOTO    208 
IF    Y(  IX 
P9=P9+1 
P(P9,0)= 
P (P9,l)= 
E=-E 
NEXT    I 
RETURN 
REP 
REM 
PRINT 
PRINT    »- 
PRINT 
FOR    1=0 
IF     1/5=1 
PRINT    P( 
GOTO    255 
PRINT    P( 
PRINT 
NEXT    I 
PRINT 
RETURN 
REM 
REM 
Z8  =  0 
E  =  C 
Kl=-1 
FOP    1=0 
IF    E=l    G 
IF    A6S(P 
ABS( 
IF    ABS(P 
E=l 

IF    ABS(P 
ABS(P( 
IF    ABS(P 
E  =  C 

K1=K1+1 
M(K1,0)= 
M  (  K 1 , 1 )  = 
M(K1,2)= 
NEXT  K 
IF  K1>=0 
PRINT  "C 
INPUT  LI 
Z8=l 
RETURN 
REM     D 
M8=M(0,1 
I9  =  C 

FOR  1=1, 
IF  ABS(M 
IF  ABS(M 
NEXT  I 
R(0,0)=M 
R(C,1)=M 


Y(0)  THEN  E=-l 

TO  LI 

)  GOTO  206C 

Y(  1-1)  GOTO  2085 

I-1  +  C8 

Y(I-l) 


Y< 1-1)  GOTO  2085 

I-1+D8 

Y(  I-1J 


***PRINT    INTERMEDIATE    PEAK    VALUES*** 

-PEAK   VALUES — « 

TO    P9 

NTd/5)     GOTC    2540 

Itl)  t 

0 

Itl) 


***DETERMINE    TWO    PEAK    VALUES*** 


TO    P9-3 

OTO    3040 

( 1,1 ) )>ABS(P( 1+1,1))     OR    ABS(P(I,1))> 

P( 1+2,1) )     GOTO    3075 

(I,l))>    ABS(P( 1+3,1) )    GOTC    3075 

(1,1 ))<ABS(P( 1+1,1)    OR    ABSfPUtUK 

1+2,1) )    GOTO    3075 

(  1,1  ))<ABS(P( 1+3,1))     GOTC    2C75 


P(I  ,0)/Sl 
P(I,1) 

I 

GOTC   3081 

HANGE    LAG    —    MAXIMUM    LAG=" 


ETERMINE    TWO    HIGHEST    PEAKS 
) 

Kl 

(  I  ,1  )  )>A6S(M8)     THEN    19=1 

<!,!))>    A8SIM8)     THEN    M8=M(I,i) 


R(C,2)= 
I8=C 
IF    19=0 
M£  =  MU6, 
IF    Kl=i 
FOR    1=18 
IF    1=19 


( 19,0) 
(19,2) 

THEN    18=1 

1) 

GOTO    3165 

+1    TC    Kl 

GOTO    3160 
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3150    IF    ABS(M(I,m>ABS  (M8)     THEN    18=1 

3155    IF    A6S(M( I, 1) )>ASS(M8)     THEN    M8=M(I,1) 

316C    NEXT    I 

3165    R(1,0)=N8 

217C    R(I,1)=M( 18, C) 

3175    R(1,2)=M(I8,23 

318C    RETURN 

3185    REM 

400C    REM  ***PRINT    TIME    DIFFERENCES*** 

4005    PRINT 

401C    PRINT 

4015    PRINT    "TIME    DIFFERENCE    BETWEEN    DIRECT    AND    SURFACE 

REFLECTED    PATHS" 
4020    PRINT    " 

4025    PRINT    TA8(9);"PEAK    VALUE" 

4027    A3  =  R(0t1  )*10C0 

403C  PRINT  USING  "       25). 25)232  MSEC", A3 

4035  PRINT 

404C  PRINT 

4045  PRINT  "TIME  DIFFERENCE  BETWEEN  DIRECT  AND  BCTTCM 

REFLECTED  PATHS" 
4050  PRINT  " 


4055  PRINT  TAE(9);"PEAK  VALUE" 

4057  A3=R(1,1 )*10C0 

406C  PRINT  USING  »       23.52222  MSEC", A3 

4065  PRINT 

407C  PRINT 

4075  RETURN 

4080  REN 

450C  REM  ***CALCULATE  SOURCE  DEPTH  AND  S-P  RANGE*** 

4505  FGR  1=1  TO  1 

451C  GC=R(0,I )*(C*R(1,  I  )  )  2 

4515  G1=C  2*R<0, I )*R( 1 , I)+4*H  2~4*H*R 

452C  G2=R<0, I )*G1 

45  2  5  G3=4*(R(0,I)=MR-H)-R*R(1,I)J 

4530  S(  I-1J=(G0-G2)/G3 

453  5  C(I-1)=2*R*S(I-1)/ (C*R (0, I ) )-C*R ( C , I)/2 

454C  NEXT  I 

4545  RETURN 

4550  REN 

500C  REM  ***PRINT  DEPTH  AND  DISTANCE  VALUES*** 

5005  PRINT  "RANGE  OF  SOURCE  FROM  RECEIVER  AND  SOURCE  DEPTH 

IN  METERS:" 
5010  PRINT  " 

— ii 

5015  PRINT  TAB(9J  ;"PEAK  VALUES" 

5020  PRINT  TAB15)  ;"DEPTH";TAB( 19)  ;"RANGE" 

503C  PRINT  USING  "    22.22222    ",S(0),D(0) 

5035  PRINT 

504C  PRINT 

5045  PRINT 

505C  RETURN 

5055  REM 

5500  REM  ***CALCULATE  DISTANCES  CF  REFLECTED  PATHS** 

5505  REM     CALCULATE  TRANSVERSE  SOURCE-RECEIVER  CISTANCE 

5510  T=SQR(D(0)  2-(R-S(C))  2) 

5515  REM     CALCULATE  DISTANCE  OF  SURFACE  REFLECTION  PATH 

552C  >=S(0J*T/(R+S(O) ) 

5525  Y=F*T/ (R+S(0)) 

5530  U=SQR(X  2+S(C)  2) 

5535  W=SQR(Y  2+R  2) 

5540  D0=U+W 

5545    REM  CALCULATE    DISTANCE    OF    BOTTCN    REFLECTICN    PATH 

5550    A=(H-S(0))*T/(2*H-R-S( C) ) 

5555    B=(h-R)*T/(2*H-R-S(0) ) 

5560    E=SQR(A    2+(H-S(0J)     2) 

55o5    F=SQR(B    2+1H-R)     2) 

5570    D1=E+F 

5575    RETURN 

5580    REM 
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600C 

REM 

***PRINT  P£TH  DISTANCES*** 

6005 
6010 

PRINT 
PRINT 

"SURFACE  REFLECTION  PATH  DISTANCE  IN  METERS:" 

6015 

PRINT 

USING  "     5122.23232",  DO 

602C 

PRINT 

6025 

PRINT 

603C 
6035 

PRINT 
PRINT 

"BOTTOM  REFLECTION  PATH  DISTANCE  IN  METERS:" 

6040 

PRINT 

USING  "     2  22. 23322" »D1 

6045 

PRINT 

605C 

PRINT 

6055 

RETCR^ 

1 

606C 

REM 

6065 

ENC 
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COMPUTER  PROGRAM  DEVERB 

r  3p  «x ^fC 3£C2fC  SQCSfC  3jS  2fE  2QC3£6  SytSQCSQCS^C  3yC  3^C  3p 3p  ?JC  SfS  *-,.  ^*  ?fC 2yC 3JC «$C  -r>  ^  t*  *  *r  -*r  * ^'5* 

C  *  DEVERB  10/21/77  * 

C  *  JEANIE    SAVAGE,     PROGRAMMER  * 

C  *  SPECIFICATIONS    BY    RICK    BOSTIAN  * 

C  *  IN    THIS    PROGRAM    DEVERBERAT I GN  * 

C  *  IS    PERFORMED    IN    THE    FRE-  * 

C  *  QUENCY    DOMAIN.  * 

C  *  LAST    CORRECTION:       12/06/77  * 

c 

INTEGER    XGRID 

INTEGERS    IZ2,IZ3,IZ4,IZ5,IZ6,YES,IZ7 
DATA    YES/'Y    •/ 

DIMENSION    I  STACK (2  0) , A(1000) ,B(  10CG) , IARYt 1C00) 
DIMENSION    X(IOOO) , Y( 1000) 

DIMENSION    FINAL ( 1  COO » 2) , I BEG( 20 ) , I  END ( 2 0 ) , I  CASE ( 20 ) 
COMMON    SF,ISF,IBM,THETA,SIGMA,GAMMA,C,D,RC0EFF,N2, 
GISIZE,DB,DS,NBLK,IZ3,IZ4 
C 

C*JU  J.  ol,  j,  j*  ..  u.4,j-,i,„  v-a.  -v  j-  gu  *u  u«  */«  -j*  u«  ^u-  .  l, 
*,*   TTTTTJrTT,T,T"iA   '.*  n*  T    -T  1*  *>»  *T  1*  "i*   -p  -i* 

C  INITIALIZATION    ROUTINE 

C 

IREAD=0 

PRINT    5  00 
500  FORM  AT ( '0*  , 'DEVERB1  ) 

PRINT    510 
510         FORMAT ( 'O* , 'NUMBER    OF    POINTS    PER    SIGNAL     (eCWER    OF    2 ) • 
£•     (15) — •) 

READ    520t     N2 
520  FQRMATU5) 

PRINT    530 
530         FORMAT (■     »,'IS    THIS    THE    FIRST    SIGNAL    OF    A    NULTIPLE', 
£«     RUN?     (Y/N)-') 

READ    540, IZ2 
540  FORMAT (Al) 

PRINT    550 
550  FORMAT!'     ', 'SAMPLING    FREQUENCY     (F9.3) — ') 

READ    560,     SF 
560  F0RMAT(F9.3) 

PRINT    5  70 
570  FORMATC     «, 'DIRECT    PATH    DISTANCE     (F9.5)  —  •) 

READ    580t     D 
580  F0RMAT(F9.5) 

PRINT    590 
590  FORMATC     ', 'SURFACE    PATH    DISTANCE    TN    METERS    (F9.5) — ') 

READ    580, DS 

PRINT    600 
600  FORMATC     ', 'BOTTOM    PATH    DISTANCE    IN    METERS     (F9.5J  — ') 

READ    580,     DB 

PRINT    610 
610         FORMATC     ', 'SURFACE    REFLECTION    TIME    IN    MSEC     (F9.5)  — '; 

READ    580,     TS 

TS=TS/1000 

PRINT    620 
620         FORMATC     ', 'BOTTOM    REFLECTION    TINE    IN    MSEC    (F9.5) — ') 

READ    580,     TB 

TB=TB/1000 

PRINT    630 
630  FORMATC     ', 'BOTTOM    REFLECTION    COEFFICIENT     (F9.5) — ') 

READ    58C,     RCOEFF 

WRITE(6,711 ) 
711  FORMATC     ','RMS    WAVE    HEIGHT    (F9.5)  — ') 

READ    580,     SIGMA 

PRINT    720 
720         FORMATC     ', 'SURFACE    ANGLE     IF    INCICENCE     (IN    RADIANS)', 
£•     (F9.5) — ' ) 

READ    580,     THETA 

PRINT    730 
730  FORMATC     ', 'BOTTOM    PHASE    SHIFT     (F9.5)  — ') 

READ    580,     GAMMA 

PRINT    640 
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640  FORMAT (•     ', 'SPEED    OF    SOUND     (F9.2) —  •) 

READ    560,     C 

PRINT  650 
650  FORMAT!"     ', 'FREQUENCY    PLOT?     (Y/N)—') 

READ    540, IZ3 

PRINT  655 
655  FORMAT!'     ','TIME     PLOT    BEFORE    FFT?    (Y/N) — •) 

READ    540, IZ7 

IF( IZ7.NE.YES)     GOTO    659 

PRINT  657 
657  FORMAT!'     ', 'NUMBER    OF    POINTS    TC    BE    PLOTTED    (15) — •) 

READ    520,     INUM1 

659  PRINT    660 

660  FORMAT(»     ','TIME    PLOT    AFTER    FFT?     (Y/N) — ') 
READ    540, IZ4 

IFUZ4.NE.YES)     GOTO    672 

PRINT    657 

READ    520,INUM2 
672  PRINT    675 

675         FORMAT (•     ',«ALL    BLOCKS?     (Y/N) — •) 

READ    540,  IZ5 

IFdZ5.EQ.YES)     GOTO    1000 

KPTR=0 

PRINT  680 
680  FORMAT!  •     «, 'SPECIFIC    BLOCKS    (12)     (INPUT    9° 

£'  FINISHED)—') 
690  READ  700,ITEMP 
700  FORMAT (I  2) 

IFdTEMP.EQ.99)     GOTO    1000 

KPTR=KPTR+1 

ISTACK(  KPTR)=dEMP 

GOTO    690 
C 

C  READ    DATA^TAPE 

£  t  *  *  *r  *  -r  *  *?  *>  *  *  *  -r  * 

C 
1000       IF(IREAD.EQ.l)     GOTO    2000 

PRINT    1005 
1005       FORMAT (•     ', 'READY    TO    READ    DATA    TAPE') 

IF(IZ2.NE.YES)     GOTO    1020 

READ    (5,1010)    KEND 
101C       F0RMAK8I6) 
102C       DO    1030    1=1, N2, 8 

ITEMP=I+7 

READ(5,1010)    (  IARY(J)  ,J  =  I ,ITEMP) 
1030      CONTINUE 
C 

C  DETERMINE    DIVISION    BOUNDARIES 

C 

PRINT    1501 
1501       FORMAT («     •  , « CONT  INU I NG    WITH    CALCLLAT IONS ■  ) 
150C       DO    1510    1=1, N2 

TIME=(I-1 )/SF 

IF(TIME.LT.TS)     GOTO    1510 

ISF=I 

GOTO  1520 
151C      CONTINUE 

ISF=N2+1 

GOTO  1540 
1520      DO    1530    1=1, N2 

TIME=( 1-1 )/SF 

IF(TIME.LT.TB)     GOTO    1530 

IBM=I 

GOTO    1550 
1530      CONTINUE 
1540       IBM=N2+1 
C 

C  ***DETERMINE    NUMBER    OF    POINTS    IN    EACH    SECTOR*** 

1550       NPT1=ISF-1 
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NPT2=IBM-ISF 

IF{ IBM.LT.ISF)     NPT2=ISF-IBM 
NPT3=N2-IBM+1 

IF( IBM.LT.ISF)     NPT3=N2-ISF+1 
C 

-Y*  -^  *^   *j*  >i»>    >j»  *(•  *p.   *»  ^»   -f»    *p  ^»  ^*    *j.    -Y"  ^*»  *r*  >f*    >A 

C  DETERMINE    BLOCK    SIZE 

C->JU    V1*    ^**    *•*-    •-*•     -J-    **•■  4»     ^X»    -J-     **•    *■*»    -Jf   «Ar    -J-     -J.  -JL.    s.1*  vJLr     -i,. 
*Tf«.  ^f*   ^j»   rf>r«   *,»    *p  -^*  .^    .-(N.   ^»    ^*   >f»   *|H  ^p.    ^f*.     *■,•   ~f.    rf+  ^f*    ^p 

c 

C  ***FIND    SMALLEST    AND    MIDDLE    NUMBER    OF    POINTS    IN    SECTOR 

1700       ISMALL=MINC(NPT1 ,NPT2) 
ISMALL=MIN0(NPT3, I  SMALL) 

IF( ISMALL.EC.NPT1 )    MI DDLE=M IN0< NPT 2, NPT3 ) 
IF(  ISMALL.EC.NPT2  )    MIDDLE=M INCH  NPT1,NPT3 ) 
IF(  I  SMALL. EC. NPT 3)     MI CDLE=M INO ( NPT1 ,NPT 2 ) 
C 

C  ***DETERMINE    NUMBER    OF    POINTS    PER    BLOCK*** 

2000    IF( IABS ( MIDDLE/2-1  SMALL). GT . (MIDDLE-ISMALL) ) 
£ISIZE=ISMALL 
IF( ISIZE.EQ.ISMALL)    GOTO    2200 
DO    2010    K=l,30 

IF( ISMALL.EC.MIDDLE/K)    GOTO    2020 

IF(  ( ISMALL.EQ.MIDDLE/K) . LT.  ( I  SMALL-MI  DDL E/ ( K+l)) } 
&GCT0    2020 
2010      CONTINUE 
2020       ISIZE=MIDDLE/K 
C 

C  ***DETERMINE    NUMBER    OF    BLOCKS    PER    SECTOR*** 

2200       NBLK1=NPT1/ISIZE 
NBLK2=NPT2/ISI ZE 
NBLK3=NPT3/ISIZF 
C 

C  ***DETERMINE    NUMBER    OF    POINTS    SKIPPED    EACH    SECTOR*** 

2300       NSKI P1=NPT1-NBLK1*ISIZE 
NSKIP2=NPT2-NBLK2*ISIZE 
NSKI P3=NPT3-NBLK3*I SIZE 
C 

C  ***PRINT    SECTOR     INFORMATION*** 

WRITE(6,2405)        ISIZE 
2405       FGRMAT( '0' t' ISIZE=» tI3) 
IT1-1 
IT2  =  2 
IT3  =  3 

WRITE(6»241C)     IT1 tNBLKl tNSKlPl 
2410    F0RMAT(  «0f  ,  'SECTOR     Ml,'     CONTAINS    't^,1     BLOCKS,1, 
$2X,I3,»     POINTS    SKIPPED1) 
WRITE(6,241C)     I T2  ,  N8LK2 , NSK I P2 
WRITE(6,241C)     I T3  ,NBLK3 , NSKI P3 
C 

C  BUILD    STACK    IF    PROCESSING    ALL    BLOCKS 

C 

2500       IFUZ5.NE.YES)     GOTO    2700 

NUMBLK=N2/I SIZE 

DO    2510    1  =  1  ,NUMBLK 

ISTACKi I)  =  I 
2510      CONTINUE 

KPTR=NUMBLK 
C 

Q  *%*^* *********************  ******* 

C  DETERMINE    BEG    6    END    POINTS    6    CASE 

Q  ********************************* 

c 

2700       16=0 
C 

C  ***BLOCKS    IN    1ST    SECTOR*** 

ITEMP=MSKIP1+1 
DO    2710    I=1,NBLK1 
16=16+1 

IBEG( 16  )=ITEMP 
IEND( I6)=ITEMP+ISIZE-1 


36 


ICASE(I6)=1 

ITEMP=IEND( I6)+l 
2710      CONTINUE 
C 
C  ***BLOCKS     IN    2ND    SECTOR*** 

ITEMP=ISF 

ITEMP1=NSKIP2 

DO  2720  I=1,NBLK2 

16=16+1 

IBEG(  16  )=ITEMP 

I END( 16 )=ITEMP+ISIZE-1 

ITEMP=I END( I6)+l 

ICASE{ I6)=2 

IF(  ISF.GT.IBMJ  ICASE(I6)=3 

IF(  ITEMP2.EQ.0)  GOTO  2720 

ITEMP=I TEMP+1 

ITEMP1=ITEMP1-1 
272C   CONTINUE 
C 
C      ***BLOCKS  IN  3RD  SECTOR*** 

DO  2730  I=1TNBLK3 

16=16+1 

IBEG( 16 )=IBN  +  (  1-1  )*ISIZE 

IF(  ISF.GT.IBM)  IBEG( I6)  =  ISF  +  ( I-1)*ISIZE 

IEND( I6)=IBEG( I6)+ISIZE-1 

ICASE( I6)=4 
2730   CONTINUE 
C 

C  PRINT  TIME  PLOT  BEFORE  FFT 

C 

2900       IF(  IZ7.NE.YES)     GOTO    3000 

IFUNUM1.GT.1000)    GOTO    2910 

NUM  =  INUM1 

KPL0TS=1 

GOTO  2920 
2910   KPLOTS=INUM1/1000+1 

NUM=1000 
2920   DC  2940  I=1,KPL0TS 

IF( ( I.EQ.KPLOTS) .AND. (I.NE. 1) )  NUM=INUM1- (KFLCTS-1 ) 
£*1000 

DC  2930  J=1,NUM 

Y( J)=FLOAT<  IARY<  I  1-1 ) *1000  + J ) ) 

X(J)=( { I-1)*1000+J-1)/SF 
2930   CONTINUE 

XGRID=NUM/1C 

IF(MOD( NUM, 10J.ME.0)  XGRI D=XGR I D+I 

CALL  PL0TDV(X,Y,NUM,XGRID,3 ,NBLK) 
2940   CONTINUE 
C 

Q  #■£■%%*%*■£■&  **■$*{■% 

C  PROCESS    BLOCKS 

c 

3000       DO    3110    1=1 ,KPTP 

NBLK=ISTACKU> 
C 
C  ***DETERMINE    BOUNDARIES    AND    CASE    CF    BLOCK*** 

IBEGO=ISEG<NBLK ) 

IENDO=I END(NBLK) 

ICASEO= ICASE(NBLK) 
C 
C      ***PROCESS  BLOCK*** 

PRINT  3015,  IBEGO,IENDO,ICASEO 
3015  FORMAT! ///•  ■ , ■  IB EG= '  ,  14, 4X , ' I END= • , 14, 4X, • ICASE= •  , 1 1 ) 

K  =  l 

DO    3020    J=IBEG0,IEND0 

A(K)  =  FLOAT(  IARY( J)  ) 

B(K)=0.0 

K  =  K+1 
3020      CONTINUE 
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302! 
3041 


305 
311 
C 

c 
c 
c 
c 

330C 


330' 


330 
33C 


331 

331 


331 
331 


3310 
3320 


CALL  SI 
I  Ft ( I  .N 

DC  3030 
SLOPE=F 
FINALU 
SLQPE=F 

F  INAL(  J 

CONTINU 
IF( (I .G 
IF(NSKI 
DO  3050 
ITEMP=I 
TEMP=FI 
FINALU 
TEMP=FI 
FINALU 
CONTINU 
CONTINU 


GNAL 

E.l) 

J  =  l 

INAL 
,1)  = 

INAL 
,2)  = 

F 

T.NB 
P2.E 
J=l 
END( 
NAL( 
TEMP 
NAL( 
TEMP 
c 


(A, 3, FINAL,IBEGO,ICASEO,ISIZE) 
.OR. (NSKIP1.EQ.0)  )    GOTC    3040 
,NSKIP1 

<IBEG(I)+1,1)-FINAL(IBEGU),1) 
-SLOPE  +  FINAL( IBEG(  I)  ,1  ) 
(IBEG(I)+1,2)-FINAL( I  BEG (I), 2) 
-SLOPE+FINAL(IBEG( I)  T2) 

LK1  +  NBLK2) .OR. (I  .LE.NBLK1  +  1  )  )     GOTO    3110 

CO)    GOTO    3110 

,NSKIP2 

NBLK1  +  D  +  1 

ITEMP+1,1  ) -FIN  ALU  TEMP- 1,1) 

,1)=FINAL( ITEMP-1,1 )+TEMP 

IT  EM  P+l«  2) -FIN  ALU  TEMP- 1,2) 

,2)=FINAL( ITEMP-1,2)+TEMP 


3330 


3340 
C 
C 
C 


IF( IZ4. 
DO  3308 
NBLK=IS 
N=ISIZE 
IBEGO=I 
IENDO=I 
K=l 

DC  3304 
A(K)=FI 
B(K)=FI 
K  =  K  +  1 
CONTINU 
CALL  CF 
DO  3306 
FINALC I 
CONTINU 
CONTINU 
IFCNSKI 
DC  3312 
FINALU 
CONTINU 
IFtNSKI 
DO  3314 
ITEMP=( 
TEMP=FI 
FINAL (I 
CCNTINU 
IF( INUM 
IFUNUM 
NUM=INU 
KPLOTS= 
GOTO  33 
KPLOTS= 
NUM=100 
DO  3340 
IF( (I  .E 
6*1000 
DC  3330 
Y(J)=FI 
Y(J)=Y( 
X(J)=( ( 
CCNTINU 
XGRID=N 
IF(M0D( 
CALL  PL 
CCNTINU 


TP  «t~  op  op  or*  op  *p  op  op  op  -op  *o~  *p  op  or*  o*  op  op  i*  or  op  *p  op  op  op 

PRINT    TIME     PLOT    AFTER    FFT 

TP  OP  Op  OP  Op  **P  Op  OP  OP  O*  OP  OP  Op  »f»  ^p  «-f»  «"■*•  *f*  <»f»  t-»  <■»»>  *■»*  IP  Op  *p 

ES)     GOTO    3500 

,KPTR 

(I) 

NBLK) 
NBLK) 

BEG0,IEND0 
Jtl  ) 
J, 2) 


NE.Y 

1=1 

TACK 

BEG( 
END( 

J=I 
NAL( 
MAL( 

E 

FT2( 
J-l 

BEGO 
E 
E 
Pl.E 

1  =  1 
,1)  = 
E 
P2.E 

1=1 
NBLK 
NAL( 
TEMP 
E 

2.GT 
2.GT 
M2 
1 

20 

INUM 
0 

1=1 
Q.KP 

J=l 
NAL( 
J)/F 

1-1) 

c 

UM/1 
NUM, 
OTDV 
E 


A,B,N,N,N,1) 
,  1ST  ZE 
+J-1,1)=A( J) 


CO)     GOTO    3313 
,NSKIP1 

<FINAL(IBEG(1),1)/(NSKIP1+1))*I 

CO)    GOTO    3315 

,NSKIP2 

1+D+l 

ITEMP+1,1 )-FINAL( ITEMP-1,1 J 

,1 )=FINAL( ITEMP-1,1 J+TEMP 

.IEND( KPTR)) INUM2=IEMD(KPTR) 
.1000)     GOTC    3310 


2/1000  +  1 

,KPLOTS 

LOTS) .AND.  (I .NE. 1) )     NUM=INUM2- ( KFLCTS-1 ) 

,  NUM 

(I-l)*1000+J,l) 
LOATUSIZE) 
*100Q+J-1)/SF 

0 

10).NE.O)     XGRID=XGRID+1 

(X,Y,NUM,XGRID,2,NBLK) 


CONCLUSION 
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[™  -r-  -r  t  t  ;p  -r  -r  "^  -■*  T> 

c 

3500       PRINT    3510 

3510       FORMATt •     •  ,'ARE    YCU    FINISHED?     (Y/N) — •) 

READ    540, IZ6 

IFUZ6.  EQ.YES)     STQP 

IREAD=1 

GCTO    672 

DE3UG    SUBCHK 

END 
C 

c 

C  ****SIGNAL    SUBROUTINE**** 

C 

SLBROUT  INE    SIGNAL (A, B , FINAL » IBEGO, ICASEO,N J 

DIMENSION    FINAL( 1000, 2) 

DIMENSION    X(IOOO) ,Y( 1000) ,Al 100  0) ,B< 1000) 

INTEGER*2    YES,IZ3,IZ4 

COMMON    SF,ISF, IBM, THETA, SIGMA, GAMMA, C,D,RC0EFF,N2, 
&ISIZE,DB,DS  ,NBLK,IZ3,IZ4 

INTEGER    XGRID 

REAL    KO 

DATA    YES/'Y    »/ 
C 

C  PERFORM    FFT 

c 

100         CALL    CFFT2( A,B,N,N,N,-1) 
C 

~r   nr  -i*  *r  -i*  -v  t  *r»  *r  n*  *i*  *r-  -r-r*  -r  -**  -v*  ^*  **» 

C  PERFORM    CORRECTIONS 

C 

500         DO    520    1=1, N 

FINAL(I8EG0+I-1,1 )=A(I) 

FINALU  BEGO  +  I-1,2  )=B(  I) 

FREQ=(I-1 )*SF/N 

PI=3. 14159 

K0=2.*P I*FREQ/C 

IF(  ICASEO.EQ.l  )     GOTO    520 

IF(  ICASE0.EC.3)     GOTO    510 
C 
C  ***CORRECTIGN    FOR    SURFACE    REFLECTION*** 

G=( (4.* PI* SIGMA* FPEQ/C)*COS (THETA ) )**2 

ITEMP={ IBEGC+I-1)-(ISF-1) 

TMAG=SQRT(FINAL( I  TEMP ,1 )**2+FINAL ( ITEMP,2)**2) 

TPHASE=ATAN2(FINAL( ITEMP, 2)  f FIN AL (  ITEMP , 1) ) 

S=D/DS*EXP(-G/2. ) *TMAG 

SPHASE=TPHASE-(DS-D)*KO-PI 

FINAL(  IBEGO  +  I-1 ,1 )=FINAL( I BEGO+ I- 1 , 1 )-S *COS (SPHASE) 

FINAL ( IBEGO  +  1-1, 2)=FINAL(  IBEGO  +  I- 1 , 2 ) -S*S IN ( SPHAS E ) 

IF(  ICASE0.EC.2)    GOTO    520 
C 

C  ***CORRECTION    FOR    BOTTOM    REFLECTION*** 

510       ITEMP=( IBEGC+I-1 )-(IBM-l) 

TMAG=SQRT  (FINAL < I  TEMP, 1 ) **2 +FIN AL ( ITEMP  ,  2)**2) 

TPHASE=ATAN2(FINAL(ITEMP,2) , F IN AL ( ITEMP , 1 ) ) 

S=RCOEFF*D/CB*TMAG 

SPHASE=TPHASE-(DB-D)*KO+GAMMA 

FINAL(IBEG0+I-1,1)=FINAL( IBEGO+I- 1,1 )-S*C CS (SPHASE ) 

FINAL( IBEGO  +  I-1, 2 )=F I NAL( IBEGO+I- 1, 2  J -S*S IN (SPHASE ) 
520         CONTINUE 
C 

C  **************************** 

C  PRINT    FREQUENCY    AND    TIME    PLOTS 

£  **************************** 

C 

1000       IF(IZ3.NE.YES)     GOTO    1020 
C 

C  ***FREQUENCY   PLOT*** 

K  =  0 
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ITEMP=N/2 

DO  1010  I=1,ITEMP 

K  =  K  +  1 

Y<K)=FINAL( IBEGO+I-1, 1)**2+FINAL(  I  8EG0  + I- 1 , 2 ) **2 

Y(K)=1D*AL0G10(Y(K)*D) 

X(K)=(K-1)*SF/N 
1010   CONTINUE 

NUM  =  K 

XGRID=NUM/1C 

IF(M0D( NUM,10) .NE.O)  XGRI D=XGR I C+ 1 

CALL  PL0TDV(X,Y,NUM,XGRID,1 ,NBLK) 
1020  RETURN 

CEBUG  SUBCHK 

END 
C 
C 

c 

C  ***#pLOT  SUBROUTINE**** 

C 

SLBROUTINE    FLOTOV{X,Y,N,XGRID,M,Ne) 

INTEGER  DfXGRIOtYGRIDtAXIS 

DIMENSION  Y(IOOO) ,C7( 101) ,0(6) , X ( 1000 ) , K AX  I  S( 51) 

DATA  IDASH/1H-/ t I  STAR /I H*/,  ID0T/1H./ 

DATA  IBAR/lHI/, IPLUS/1H+/, I BLANK/1H  /,IX/1HX/ 

AXIS=51 

YGRID=6 

XGRID=XGRID+1 
2120   N1=N-1 

Y6=Y(1) 

Y1=Y(1) 

DO  2200  1=1 ,N 

I F(Y6-Y(I ) .GE.O.O)  GOTO  2180 

Y6  =  Y( I) 
2180   IF( Yl-Y( I  ).LE.O.O)  GOTO  2200 

Y1  =  Y(I  ) 
2200   CONTINUE 

S=Y1*(AXIS-1)/(Y6-Y1) 

X1=X(1) 

X10=X<N) 

0( 1)=Y1 

0(6)=Y6 

IIX=XGRID-1 

DO  2410  1  =  1, IIX 

C7( I)=X(( I-l)*10+l) 
2410   CONTINUE 

C7(XGRID)=C7(XGRID-1  H10*  (  X  {  2  )-X  (  1 )  ) 

IF(N.EQ.(XGRID-1 )*10)  C7(XGRID)=X(N) 

I  IY=YGRID-1 

DC  2440  1=2, IIY 

0(I)=(FL0AT(I-1 )*(Y6-Y1 )/FLO AT( YGR ID-1 ) )+Yl 
2440   CONTINUE 

WRITE(6,2460) 
2460   FORMAT<///,«  ») 

IF(M.NE.l)  GOTO  2485 

PRINT  2470, NB 
2470   FORMAT ( «0' ,32X, • BLOCK* , IX, I  2  ) 

2485  IF(M.EQ.l)     PRINT    2486 

2486  FCRMAT(»     •,27X,»DB,,S    VS.    FREQUENCY') 
IFCM.EQ.2)     PRINT    2488 

2488       FORMATP     ', 23X, • VOLTAGE    VS.     TIME     (AFTER    FFT)») 
IF(M.EQ.3)     PRINT     2487 

2487  FORMAT!1     •  ,  23X  ,  '  VOLTA  GE    VS.     TIME     (BEFORE    FFT)M 
WRITE(6,2500)     ( 0(  I ) , I =1 , YGR I  0 ) 

2500      F0RMAT(9X,1H1PE10.2)  ) 

S1=(X10-X1)/10.0*(XGRID-1) 

D=l 

Ll=l 

L  =  l 

I2=IFIX(-S+1.5) 

ITEMP=( XGRID-1 )*10+1 

DO    2900     I1=1,ITEMP 

IF(N.LT.ll)    GOTO    2510 
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2510 
265C 
268C 


2720 
2725 


2760 
2780 

2810 


2860 
2865 
287C 

2900 
2910 


C 
C 
C 
C 


10 


2720 

(0.0.LE.Y6) ) 


KAXIS( IZ )=ICCT 


YTEMP=( Y( II )*FLOAT(YGRID-l) *10. 0/ ( Y6-Y1 ) )-S 

IY=IFIX(YTFMP+1 .5) 

IF(Ll.GT.L)  GOTO  2760 

DO  2650  IP=1,AXIS 

KAXIS< IP)=IDASH 

CONTINUE 

DO  2680  1=1, AXIS 

KAXIS( I )=IPLU5 

CONTINUE 

IF(N.LT.Il)  GOTO 

IFUY1.LE.0.0) .AND 

KAXIS( IY)=ISTAR 

WRITE(6T2725)  C7(D), (KAXISi  J) ,J=1,AXIS) 

FCRMAT( 1PE13.2,2X,115A1) 

L1=L1+10 

D  =  D+1 

GCTO  2370 

DO  2780  IP=1,AXIS 

KAX1S( IP)=1BLANK 

CONTINUE 

DC  2810  1=1, AXIS, 10 

KAXIS( I  )=IBAR 

CONTINUE 

IF(N.LT.Il)  GOTO  2860 

IF((Y1.LE.0.0).AND.(0.0.LE.Y6))  KAXISI I  I )  =  ICOT 

KAXIS( IY)=ISTAR 

WRITE (6, 2 86 5)  (KAXIS(JJ,J  =  1,AXIS) 

FORMAT(  15X,  115A1 ) 

L  =  L+1 

IFU.GT  .(XGRID-l)*10  +  2)    GOTO    291C 

CONTINUE 

RETURN 

CE8UG  SUBCHK 

END 

***FOURIER    TRANSFORM    SUBROUTINE*** 
SUBROUTINE    CFFT 2 < A , B , NTOT , N , NSP AN  ,  I SN ) 
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COMPUTER  PROGRAM  TDEVERB 

Q  -if.  i}  *  # *  * ^  ^0?  ^r^*;^  ^^^^jji^i:*^^^^-.^^^*^^:^^** 

C  *  TDEVERB  12/07/77  * 

C  *  JEANIE    SAVAGE,     PROGRAMMER  * 

C  *  SPECIFICATIONS    BY    RICK    BOSTIAN  * 

C  *  IN    THIS    PROGRAM    DEVERBERAT I  ON  * 

C  *  IS    PERFORMED    IN    THE    TIME  * 

C  *  DOMAIN.  * 

C  *  LAST    CORRECTION:       12/08/77  * 

c 

INTEGER    XGRID 

INTEGERS     IZ2,IZ3,IZ4,IZ5,IZ6,YES,IZ7 

DATA    YES/'Y    '/ 

DIMENSION    I  STACK (20) , A( 1000 ) , B( 1000) , IARY( 1  COO) 

DIMENSION    X(IOOO) ,Y(1000) 

DIMENSION    FINAL (1000)  ,  IBEG<  20  J ,  I  END (20) , ICASE(20) 
C 

C  INITIALIZATION    ROUTINE 

CJ »  «J-  JL»  »t»  «Jf  •*!»  *JU  -J-.  »U  0»  J*  -X  «J»  O-  -JL.  «JL  >*-  ^  -  -Jt*-  -A»  -J-  >U  «J  , 
■rr*  -n*  o*  •^r*-r-  *r  Jj*  -i*  «s*  ^-  -r  ^s-^^^s  •rxT'r^'r  ^  -r 

C 

IREAD=0 

PRINT    500 
500  FORMAT! '0' ,'DEVERB' ) 

PRINT    510 
510  F0RMAT( '0' , 'NUMBER    OF    POINTS    PER    SIGNAL     (PCWER    0C    2)\ 

£•       (15) — •)  ' 

READ    520,     N? 
520  F0RMAT(I5) 

540         FORMAT(Al) 

PRINT    550 
550         FORMATS     ', 'SAMPLING    FREQUENCY     (F9.3) —  ') 

READ    560t     SF 
560  FORMAT(F9.3) 

PRINT    5  70 
570  FORMAT (■     ', 'DIRECT    PATH    DISTANCE     (F9.5) — •) 

READ    580,     0 
580         F0RMAT(F9.5) 

PRINT    5  90  . 

590         FORMAT  (■     ', 'SURFACE    PATH    DISTANCE    IN    METERS    (F9.5)  —  'J 

READ    580, DS 

PRINT    600 
600  FORMAT (•     ','BOTTCM    PATH    DISTANCE     IN    METERS     (F9.5) — ») 

READ    580,     DB 

PRINT    610  .     , 

610         FORMAT  (•     ', 'SURFACE    REFLECTION    TIME     IN    MSEC    (F9.5)  — ■) 

READ    580,     TS 

TS=TS/10O0 

PRINT    620 
620  FORMAT ('     ','BOTTCM    REFLECTION    TINE    IN    MSEC     (F9.5) — ') 

READ    580,     TB 

TB=TB/1000 

PRINT    630 
630         FORMAT (•     ','BOTTCM    REFLECTION    CCEFFICIENT     (F9.5) — • ) 

READ    580,     RCOEFF 

WRITE(6  ,711 ) 
711  FORMAT(«     •  ,'RMS    kAVE    HEIGHT    (F9.5)  — ') 

READ    580,     SIGMA 

PRINT    720 
720         FORMAT (•     ', 'SURFACE    ANGLE    OF    INCIDENCE     (IN    RADIANS)1, 
£•        (F9.5)— •) 

READ    580,     THETA 

PRINT    730 
730         FORMAT (•     «, 'BOTTOM    PHASE    SHIFT     (F9.5)  — ') 

READ    530,     GAMMA 

PRINT    640 
640         FORMAT («     ', "SPEED    OF    SOUND     (F9.3)  — ') 

READ    560,     C 

PRINT    645 
645    FCRMATt «     ', 'BLOCK    SIZE     (15) — ') 

READ    520, ISIZE 

PRINT    6  50 
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c 
c 
c 

c 
c 


650         FORMAT!'     ', 'FREQUENCY    PLOT    BEFORE    CORRECT  I CNS?f T 
£•        (Y/N ) — ■  ) 

REAO    540,  IZ3 

PRINT    652 
652    FCRMAT(  •     ', 'FREQUENCY    PLOT    AFTER    CORRECTIONS?    (Y/N) — ') 

READ    540,  IZ2 

PRINT    655 
655  FORMAT {»     ','TIME    PLOT    BEFORE    CCRRECTIONS?     (Y/N) — ') 

READ    540,  IZ7 

IF(  IZ7.NE.YES)     GOTO    659 

PRINT    657 
657         FORMAT!'     », 'NUMBER    OF    POINTS    TO    BE    PLOTTED    (15) — •) 

READ    520,     INUM1 

659  PRINT    660 

660  FORMAT! •     ','TIME    PLOT    AFTER    CORRECTIONS?     (Y/N) — ') 
READ    540, IZ4 

IFdZ4.NE.YES)     GOTO    672 

PRINT    657 

READ    520, INUM2 
672         PRINT    675 
675         FORMAT (•     ','ALL    BLOCKS?     (Y/N) — •) 

READ    540, IZ5 

IF(  IZ5.EQ.YES)     GOTO    1000 

KPTR=0 

PRINT    680 
680         FORMAT!'     »,  'SPECIFIC    BLOCKS    (12)     (INPUT    99    WFEN    FIMS 
690  READ    700,ITEMP 

700  F0RMATU2) 

IF(  ITEMP.EQ.99  )     GOTO    1000 

KPTR=KPTR+1 

ISTACK(KPTR)=ITEMP 

GOTO    690 

3|c  jjs  ^f  ^k  ^  ^  -*  ^  -/  ^~  ^' ;  ?£  ^c  ^Jt  ?"t 

READ    DATA    TAPE 

■)»■     mt  %i*  «j*^  iJU  *JL>  tjm  O*  *A»    <-l-  -JL»   0»  *JL>  -J* 

*y        J-|k    »^»     ^p    >p     ,>p      ^,41    *rv    -(-.      ^p     ^f.      *[-     *ya    -.» 


1000 

1005 
1010 
1020 


C 
C 

c 

c 

c 


1030 


1501 
1500 


1510 
1520 


C 

C 

c 
c 
c 


1520 
1540 


IF( IREAD.EQ.l )     GOTO    2500 
PRINT    1005 

FORMAT (•     ', 'READY    TO    READ    DATA    TAPE') 

F0RMAT(8I6) 

DO    1030    1  =  1, N2, 8 
ITEMP-I+7 
READ(5,1010)    (  IARY( J)  ,J  =  1 ,ITEMP) 

CONTINUE 

£  *  *  **  ***  *  *  *  ***  ***  **  **  *  *  *  *  *  *  *  * 

DETERMINE    DIVISION    BOUNDARIES 

*****£************* ********3* 

PRINT    1501 

F0RMAT(«     ', 'CONTINUING    WITH    CALCULATIONS') 

DO    1510    1=1, N2 
TIME=(I-1 )/SF 
IF(TIME.LT.TS)     GOTO    1510 
ISF=I 
GOTO    1520 

CONTINUE 
ISF=N2+1 
GOTO    1540 

CO    1530    1=1, N2 
TIME=(I-1 )/SF 
IF(TIME.LT.TB)     GOTO    1530 
IBM=I 
GOTO    2500 

CONTINUE 

IBM=N2+1 


*  **  **  ****  *  *  ********  **  *  ***  **  4  *****  ** 

BUILD    STACK     IF    PROCESSING    ALL    BLOCKS 

************************************ 


43 


c 

250C       IF( IZ5.NE.YES)     GOTO    2900 
NUMBLK=N2/ISIZE 

DO    2510    I=1,NUMBLK 

ISTACK( I)  =  I 
251C      CONTINUE 

KPTR=NUMBLK 
C 

'("■  t  ;r  -v-  - >  <^  *i*  t»-v  Is  -r»  -v  ^  "r  *r  *r  f  *t*  -r  *r  -i*  n»  -t»  -i-  -nr*  -r  t  ^r  -r  -*r  -t*  t»  -i*  -v* 

C  PRINT    TIME     PLOT    BEFORE    CORRECTIONS 

t*  -T  fl*  •^*tt  -r  i*  n*  -T  i*  T*  J?  >-  -t*  ^  *F  T*  *P  -T-  -r»'(*  *r  "P  -r*  *r  ^*-r  'r  >?  *t*  -r-  -r  f*  -P 

C 

2900       IF< IZ7.NE.YES)     GOTO    3000 

IF( IZ5.NE.YES)     GOTO    2950 
C 
C  ***TIME     PLOT    FOP     ENTIRE    SIGNAL*** 

IF(INUMl.GT.lOOO)     GOTO    2910 

NUM=INUM1 

KPL0TS=1 

GOTO    2920 
2910       KPL0TS=INUM1/1000+1 

NUM=1000 
2920       DO    2940    I=l,KPLOTS 

IF(( I.EQ.KPLOTS) .AND. (I.NE. 1) )     NUM-INUMl-< KFLCTS-1 ) * 
&1C00 

DG    2930    J=1,NUM 

Y< J)=FL0AT{  IARY(  (  1-1)  *1000  +  J)  ) 

X(J)=( ( I-1)*1000+J-1)/SF 
2930      CONTINUE 

XGRID=NUM/10 

IF(MOD( NUM, 10).NE.G)     XGRID=XGRI C+l 

NBLK=-I 

CALL    PL0TDVIX,Y,NUMTXGRID,3,N3LK) 
2940      CONTINUE 

GCTO    3000 
C 

C  ***TIME    PLOT    FOR    INDIVIDUAL    BLOCKS*** 

2950    DC    2970    1=1 ,KPTR 

NELK=ISTACK(I ) 

DO    2960    J=1,ISIZE 

Y( J)=FLOAT( IARY(NBLK-1)*ISI ZE+J) 

X( J)=( (NBLK-1)*ISIZE+J-1 )/SF 
296C    CCNTINUE 

XGRID=ISIZE/10 

IF(MOD(  ISIZEtlOJ.NE.O)     XGRI D=XGR I C  +  l 

CALL    PLCTDViX, Y, ISIZE,XGRID,3»NBLK) 
2970       CONTINUE 
C 

c 
c 

C  PERFORM    CORRECTIONS 

C 

3000    DC    3010    1  =  1  ,N2 

FPEQ=U-1 )*SF/ISIZE 
PI=3. 14159 

G=((4.*PI*SIGMA*FREG/C)*C0S<THETA) )**2 
FINAL( I )=FLCAT(I ARY( I) ) 
ITEMP=I-ISF+1 
ITEMP2=I-I3M+1 

IF(I.GE.ISF)    FINAL (I ) =FINAL ( I ) + EXP (-G/2 . )*0*DS* 
&FINALUTEMP) 

IF( I.GE.IBM)    FINAL (I )=FINAL(I )-RCCEFF*D*DB* 
£FINAL< ITEMP2) 
301C    CCNTINUE 
C 

C  PRINT    TIME     PLOT    AFTER    CORRECTIONS 

c 

IF(  IZ4.NE.YES)     GOTO    3400 
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IF( IZ5.NE.YES)  GOTO  3350 
C 
c      **-TlME  PLOT  FOR  ENTIRE  SIGNAL*** 

IF( INUM2.GT.1000)  GOTO  3310 

NUM=INUM2 

KPL0TS=1 

GOTO  3320 
3310  KPLOTS=INUM2/1000+1 

NUM=1000 
3320  DC  3340  I=1,KPL0TS 

IF( ( I.EQ.KPLOTS) .AND.  (I.NE.l) )  NUM=INUM 2-( K PLOTS- 1 ) 
£1000 

DO  3330  J=1,NUM 

Y( J)=FI NAL( ( 1-1 )*1000+J) 

X(J)=( ( I-1)*1000+J-1)/SF 
3330    CONTINUE 

XGRIO=NUiM/lC 

IF (MQO( NUM,10).NE.O)     XGRID=XGRI C+l 

NBLK=-I 

CALL  PL0TDV(X,Y,NUM,XGRID,2 ,NBLK) 
3340  CONTINUE 

GCTO  3400 
C 

C      **^time  PLOT  FOR  INDIVIDUAL  BLOCKS*** 
3350  DO  3370  1=1 tKPTR 

NeLK=ISTACK(I ) 

DO  3360  J=1,ISIZE 

Y( J)=FINAL( <NBLK-1)*ISIZE+J) 

X( J)=( (NBIK-1)*ISIZE+J-1)/SF 
3360  CONTINUE 

XGRID=ISIZF/10 

I F(MOD(  ISIZE,10).NE.O)  XGRI D=XGRI C+l 

CALL  PLCTDV(X, Y,I SIZE,XGRID,2,NELK) 
3370  CONTINUE 
C 
C 

C  PRINT  FREQUENCY  PLOTS 

C-J-  J-  J*  J.  V-  J-  J-  -J-  -■-  -  K,   -J-  ,  >  O .  ^«,  O*  -'-  -X,  .JU  -J U  ,1. 
~r  *(-   J|S  ~f   -*.  *(+   J?  ffi  ^p.  -,*  ^  y-  *,»  ^»  -,*  *,»  -r  -,----,-  ,,. 

C 

3400  IFUZ3.NE.YES)  GOTO  3450 
C 
C      ***FREQUENCY  PLOT  BEFORE  CORRECT ICNS*** 

DO  3430  1=1 ,KPTR 

NBLK=ISTACK(I ) 

DC  3410  J=1,ISIZE 

A( J)=FLOAT< IARY( (NBLK-1)*ISIZE+J)) 

B(J)=0.0 
3410  CONTINUE 

CALL  CFFT2( A,B,ISIZE, ISIZE,ISIZE,-1) 

ITEMP=ISIZE/2 

DO    3420    J=1,ITFMP 

Y( J)=A( J)**2+B(J )**2 

Y( J)=10*AL0G10<Y( J) ) 

X( J)=( J-l )*SF/ISIZE 
3420    CONTINUE 

XGRID=ITEMP/10 

IF    (MOD( ITEMP,10) .NE.O)    XGR ID=X GR  I  C+l 

CALL    PLOTDV(X,Y,ITEMP  ,XGRID , 1 ,N ELK ) 
3430    CONTINUE 
C 

C  ***FREQUENCY    PLOT    AFTER    CORRECT  I  C!\S*** 

3450    IF(IZZ.NE.YES)     GOTO    3500 

DO    3480    I=1,KPTR 

NBLK=ISTACK(I ) 

DC  3460  J=1,ISIZE 

A(J)=FINAL((NBLK-1)*ISIZE+J) 

B{  J)  =  0.0 
346C    CONTINUE 

CALL    CFFT2( A,B, I  SIZE,  I  SIZE,  ISI 2E t-l) 

ITEMP=ISIZE/2 

DO    3470    J=1,ITEMP 
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Y( J)=A( J)**2+B< J)**? 
Y( J)=10.*ALCG10(Y( J) ) 
X( J)=(J-1)*SF/ISI ZF 
347C  CCNTINUE 

XGRID=I TEMP/10 

IF(MOD(  IT  EM P, 10)  .NE.O)  XGRI D=XGR  I  C  +  l 

CALL    PLOTDV(X,Y,ITEMP  ,XGRID , 4 , N ELK ) 


348C 
C 
C 

CCNTINUE 

*i*  3{e  ^r  *r  ^  3?«.*~  !?(  >p 

C 

CONCLUSION 

C 

c 

3500 

^<3j:^t^t^;^!^£^rf:^t 

PRINT  3510 

3510 

FORMAT* •  ','ARE  YOU  FINISHEO? 

READ  540, IZ6 

IF( IZ6.EQ.YES)  STOP 

IREAD=1 

GCTO  672 

DEBUG  SUBCHK 

END 

C 
C 
C 
C 
C 
C 


2120 


2180 
2200 


(Y/N)--' ) 


2410 


2440 
2460 

2465 

2466 
2470 
2485 
2486 

2438 


SLBROUTIN 

INTEGER  D 

DIMENSION 

DATA  IDAS 

DATA  IBAR 

AXIS=51 

YGRID=6 

XGRID=XGR 

N1=N-1 

Y6=Y(1) 

Y1=Y(1 ) 

DO  2200  I 

IF(Y6-YU 

Y6=Y(I) 

IF(Y1-Y(I 

Y1=Y( I ) 

CONTINUE 

S=Y1*( AXI 

X1=X(1) 

X10=X(N) 

0(1)=Y1 

G(6)=Y6 

IIX=XGRID 

DO  2410  I 

C7(I)=X( ( 

CONTINUE 

C7(XGRID) 

IF(N.EQ.( 

IIY=YGRID 

DO  2440  I 

0(I)=(FLO 

CONTINUE 

WRITE( 6,2 

FORMATt // 

IFINB.GT. 

NE=-NB 

PRINT  2  46 

FORMAT( '0 

GOTO  2485 

PRINT  247 

FORMAT! »0 

IFiM.EQ.l 

FCRMAT( « 

I F(M.EQ.2 

FORMATt  * 


*#**Pj_OT  SUBROUTINE**** 

E  PLOTDV(XtY,NtXGRID,M,Ne) 
,XGRID,YGRIDTAXIS 

Y(IOOO) ,C7( 101) ,0(6) tX(lOOO) ,KAXIS(51) 
H/1H-/, I STAR/1H*/,  I  DOT /IF./ 
/1HI/, IPLUS/1H+/,! BLANK/ 1H  /,IX/1H>/ 


ID+1 


=  1  ,N 

).GE.O.O)  GOTO  2180 

J.LE.O.O)  GOTO  2200 


S-1)/(Y6-Y1) 


-1 
■l.IIX 

I-l)*10+l) 

=C7(XGRID-1 ) +10*(X(2)-X(1) ) 
XGRID-1)*10)  C7<XGRID)=X(N) 
-1 

=  2tHY 

AT(I-1 )*(Y6-Y1)/FL0AT(YGRID-1) ) +Y1 

46C) 
/»  '  '  ) 

0)  GOTO  2466 

5f  NB 

' t  32X, • PLOT' ,1X, 12) 

0,NB 

• ,32X, ' BLOCK' ,1X,I2) 

)  PRINT  2436 

■,17X,'DB"S  VS.  FREQUENCY  (BEFORE  CORRECTICN 

)  PRINT  7488 

•t20Xt 'VOLTAGE  VS.  TIME  (AFTER  CORRECT  I ONS )•  ) 
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IF(M.EQ.3)     PRINT    2487 
2487      FORMAT!1     • , 20X, ' VOLTA GE    VS.     TIME     (BEFORE    CORRECTIONS ) ■) 

IF1M.EQ.4)     PRINT    2489 
2489    FCRMAT(«     ',17X,'DE"S    VS.    FREQUENCY    (AFTER     ', 
^'CORRECTIONS)' ) 

WRITE(6,2500)    ( 0(  I  ) ,  I =1 , YGR I D ) 
2500      FORiMAT(9X,ll(lPE10.2)  ) 

S1=(X10-X1)/10.0*(XGRID-1) 

D-l 

Ll  =  l 

L  =  l 

IZ-IFIX(-S+1.5) 

ITEMP=( XGRIC-1 )*1C+1 

DO    2900    I1=1,ITEMP 

IF(N.LT.Il)    GOTO    2510 

YTEMP=( Y( II  )*FL0AT(YGRID-1)*10.0/(Y6-Y1)  )-S 

IY=IFIX(YTEMP+1.5) 
2510  IF(Ll.GT.L)  GOTO  2760 

DO  2650  IP=1,AXIS 

KAXIS( IP)=IDASH 
2650   CONTINUE 

DC  2680  1=1, AXIS, 10 

KAXISU  )  =  IPLUS 
2680   CONTINUE 

IF(N.LT.Il)  GOTO  2720 

IF((Y1.LE.0.0).AND.(0.0.LE.Y6))  KAXISU  Z)=ICOT 

KAXIS( IY)=ISTAR 
2720   WRITE(6,2725)  C 7 ( D ) , ( KAX I S ( J } , J  =  1 , AX  IS ) 
2725   FCRMATt  1PE  12.2 , 2X , 11 5A1 ) 

L1=L1+10 

D  =  D+1 

GCTO  2370 
2760   DC  2780  IP=1,AXIS 

KflXIS( IP)=IBLANK 
2780       CONTINUE 

DC    2810    1  =  1, AXIS, 10 

KAXIS(  I  )  =  IBAR 
2810   CONTINUE 

IF(N.LT.Il)  GOTO  2860 

IF((Y1.LE.0.0).AND.(0.0.LE.Y6) )  KAXI S ( I Z )  =  1  COT 

KAXIS( IY)=I STAR 
286C   WRITE(6,2865)  ( KA XI S ( J) , J  =  l , AX  I S J 
2865   FORMAT!  15X, 115A1 ) 
2870   L=L+1 

IF(L.GT.(XGRID-1 )*10+2)  GOTO  2910 
2900   CONTINUE 
291C   RETURN 

CEBUG  SUBCHK 

END 
C 
C 
C  ***FOUPIER  TRANSFORM  SUBROUTINE*** 


C 


SUBROUTINE  CFFT2( A, B , NTOT ,N , NSP AN , ISN) 
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ing spectral  source 
levels  of  marine  mam- 
mals in  coastal  waters. 
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